
#include <navierstokes.hpp>
#include <feel/feelfilters/geotool.hpp>
#include <iostream>

//#include <feel/feelfilters/gmshhypercubedomain.hpp>
#include <feel/feelfilters/gmshsimplexdomain.hpp>

#include <research/life-workspace/Fluid-Structure/functionSup.cpp>


namespace Feel {

namespace FSI {

template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
NavierStokes<Dim,Order,Entity>::NavierStokes( application_ptrtype __appli )
            :
            M_lifeApplication(__appli),
            M_backend( backend_type::build( this->application()->vm() ) ),
            meshSize( this->application()->vm()["fluide-hsize"].template as<double>() ),
            mesh_fluid( this->createFluidMesh( meshSize,"domainFluid" )),
            mesh_reference_fluid( this->createFluidMesh( meshSize,"domainFluidRef" )),
            mesh_struct_top(this->createStructMesh( meshSize,1.0, "boundaryMoveH" ) ),
            mesh_struct_bottom(this->createStructMesh( meshSize,0.0, "boundaryMoveB" ) ),
            Sh_top(new structure_space_type( mesh_struct_top )),
            Sh_bottom(new structure_space_type( mesh_struct_bottom )),
            // reference boundary description
            M_bc_reference_top(new element_structure_type( Sh_top, "bc_reference_top" ) ),
            M_bc_reference_bottom(new element_structure_type( Sh_bottom, "bc_reference_bottom" ) ),

            M_Xh(new space_fluid_type( mesh_fluid ) ),
            M_Xh_ref(new space_fluid_type( mesh_reference_fluid ) ),
            U_fluid( new element_fluid_type(M_Xh,"U")),
            M_Ufluid_ref( new element_fluid_type(M_Xh_ref,"Uref")),
            M_identity_ale(new element_fluid_velocity_type( M_Xh->template functionSpace<0>(), "identity_ale" ) ),
            M_mesh_velocity(new element_fluid_velocity_type( M_Xh->template functionSpace<0>(), "mesh_velocity" ) ),
            M_bdf_fluid ( new bdf_type( this->application()->vm(),M_Xh,"fluid" )),
            M_depl_top(new element_structure_type( Sh_top, "bc_top" )),
            M_depl_bottom(new element_structure_type( Sh_bottom, "bc_bottom" )),
            M_aleFactory( new ALE< convex_type >( std::make_pair(0,4), mesh_reference_fluid, this->application()->vm() ) ),
            M_aux(new ale_map_element_type( M_aleFactory->getDisplacement().functionSpace(), "new ale map" ) ),
            M_bdf_ale_map ( new bdf_ale_map_type(this->application()->vm(),
                                                 M_aleFactory->getDisplacement().functionSpace(),
                                                 "ale_map") ),
            M_bdf_ale_identity ( new bdf_ale_velocity_type( this->application()->vm(),
                                                            M_Xh->template functionSpace<0>(),
                                                            "velocity fluid" ) ),
            M_exporter( Exporter<mesh_type>::New( this->application()->vm(), this->application()->about().appName() ) ),
            M_exporter_ref( Exporter<mesh_type>::New( this->application()->vm(), this->application()->about().appName()+"_ref" ) )
         {
             M_time = __appli->vm()["bdf-time-initial"].template as<double>();
             M_cpt_export=0;
            //Masse volumique
             rho = this->application()->vm()["rho"].template as<double>();
             M_mu = this->application()->vm()["mu"].template as<value_type>();
             penalbc = this->application()->vm()["bccoeff"].template as<value_type>();
             Vitesse_Entree = this->application()->vm()["entree-v"].template as<value_type>();
             Pression_Sortie = this->application()->vm()["sortie-p"].template as<value_type>();
             M_weakCL = this->application()->vm()["weakCL"].template as<bool>();
             M_isStat = this->application()->vm()["stationnaire"].template as<bool>();
             M_isMoveDomain = this->application()->vm()["movedomain"].template as<bool>();
             M_isFixeBoundary = false;

             *M_bc_reference_top = vf::project( Sh_top, elements(Sh_top->mesh()), vf::constant(1.0) );
             *M_bc_reference_bottom = vf::project( Sh_bottom, elements(Sh_bottom->mesh()), vf::constant(0.0) );
             M_referencePolyBoundarySet.push_back(*M_bc_reference_bottom);
             M_referencePolyBoundarySet.push_back(*M_bc_reference_top);


             M_LookComment = this->application()->vm()["fluide-verbose"].template as<bool>();//false;
#if 0
            // set up the non linear solver
            M_backend->nlSolver()->residual = boost::bind( &NavierStokes<Dim, Order, Entity>::updateResidual,
                                                           boost::ref( *this ), _1, _2 );
            M_backend->nlSolver()->jacobian = boost::bind( &NavierStokes<Dim, Order, Entity>::updateJacobian,
                                                           boost::ref( *this ), _1, _2 );

#endif
            std::cout << "constructeur NS finsih\n";

         }


template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
typename NavierStokes<Dim,Order,Entity>::mesh_ptrtype
NavierStokes<Dim,Order,Entity>::createFluidMesh( double meshSize,std::string name )
{
    if (M_LookComment) std::cout << "[NavierStokes] : createFluidMesh start\n";


    this->application()->changeRepository( boost::format( "%1%/%2%/P%3%/h_%4%/" )
                                           % this->application()->about().appName()
                                           % convex_type::name()
                                           % Order
                                           % this->application()->vm()["fluide-hsize"].template as<double>()
                                           );


    //Geometry
    GeoTool::Node x1(0,0);
    GeoTool::Node x2(4,1);
    GeoTool::Rectangle R( meshSize,"LEFT",x1,x2);
    R.setMarker(_type="line",_name="Entree",_marker4=true);
    //R.setMarker(_type="line",_name="Paroi",_marker1=true,_marker3=true);
    R.setMarker(_type="line",_name="ParoiB",_marker1=true);
    R.setMarker(_type="line",_name="ParoiH",_marker3=true);
    R.setMarker(_type="line",_name="Sortie",_marker2=true);
    R.setMarker(_type="surface",_name="OmegaFluide",_markerAll=true);

    GeoTool::Node x_centre(1,0.6);
    GeoTool::Node x_bord(1.1,0.6);
    GeoTool::Circle C( meshSize/3.,"LEFT2",x_centre,x_bord);
    C.setMarker(_type="line",_name="Paroi",_markerAll=true);
    C.setMarker(_type="surface",_name="OmegaFluide",_markerAll=true);

    GeoTool::Node x_centre2(2.5,0.3);
    GeoTool::Node x_bord2(2.7,0.3);
    GeoTool::Circle C2( meshSize/3.,"disque2",x_centre2,x_bord2);
    C2.setMarker(_type="line",_name="Paroi",_markerAll=true);
    C2.setMarker(_type="surface",_name="OmegaFluide",_markerAll=true);


    if (M_LookComment) std::cout << "[NavierStokes] : createFluidMesh finish\n";

    return R.createMesh<mesh_type>(name);



} // NavierStokesStokes::createFluidMesh




template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
typename NavierStokes<Dim,Order,Entity>::mesh_structure_ptrtype
NavierStokes<Dim,Order,Entity>::createStructMesh( double meshSize,double __ypos, std::string name )
{
    if (M_LookComment) std::cout << "[NavierStokes] : createStructMesh start\n";

    this->application()->changeRepository( boost::format( "%1%/%2%/P%3%/h_%4%/" )
                                           % this->application()->about().appName()
                                           % convex_type::name()
                                           % Order
                                           % this->application()->vm()["fluide-hsize"].template as<double>()
                                           );

#if 0
    mesh_structure_ptrtype struct_mesh( new mesh_structure_type );
    //GmshHypercubeDomain td(1,1,1,true);
    GmshSimplexDomain td(1,1);

    td.setCharacteristicLength( meshSize );
    td.setX( std::make_pair( 0, 4 ) );
    std::string fname = td.generate( "boundary_description" );
    ImporterGmsh<mesh_structure_type> struct_import( fname );
    struct_mesh->accept( struct_import );

    if (true) std::cout << "[NavierStokes] : createStructMesh finish\n";

    return struct_mesh;
#endif

    GeoTool::Node x1(0,__ypos);
    GeoTool::Node x2(4,__ypos);
    GeoTool::Line R( meshSize,"kokok",x1,x2);

    if (M_LookComment) std::cout << "[NavierStokes] : createStructMesh finish\n";

    return R.createMesh<mesh_structure_type>(name);
}



template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
void
NavierStokes<Dim, Order, Entity>::solve( sparse_matrix_ptrtype const& D,
                                   element_fluid_type& u,
                                   vector_ptrtype const& F,
                                   bool is_sym )
{
    vector_ptrtype U( M_backend->newVector( u.functionSpace() ) );
    M_backend->solve( D, D, U, F, false );
    u = *U;
} // NavierStokes::solve

template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
void
NavierStokes<Dim, Order, Entity>::exportResults( element_fluid_type& U )
{
    M_exporter->step( 0 )->setMesh( U.functionSpace()->mesh() );
    M_exporter->step( 0 )->add( "u", U.template element<0>() );
    M_exporter->step( 0 )->add( "p", U.template element<1>() );
    M_exporter->save();
} // NavierStokes::export


template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
void
NavierStokes<Dim, Order, Entity>::exportResultsInstat( mesh_ptrtype __mesh,element_fluid_type& U,element_fluid_velocity_type & mesh_velocity,
						       ale_map_element_type & aux,double time )
{

    if (true) std::cout << "[NavierStokes] : Export NS start\n";


    //this->application()->changeRepository( boost::format( "%1%/%2%/P%3%/h_%4%/Export/" )
    this->application()->changeRepository( boost::format( "%1%/%2%/P%3%/h_%4%/" )
                                           % this->application()->about().appName()
                                           % convex_type::name()
                                           % Order
                                           % this->application()->vm()["fluide-hsize"].template as<double>()
                                           );

#if 0
    //getSolution()->zero();

    //element_fluid_ptrtype Stress(new element_fluid_type (M_Xh,"Stress"));
    
    //auto Stress = getSolution();
    //element_fluid_0_type uStress = Stress->template element<0>();

    //modifVec222(uStress,Stress, -vf::oneY()/*Sigmav*N()*/ , "ParoiB" );
    //modifVec222(getSolution()->template element<0>(),getSolution(), -vf::oneY()/*Sigmav*N()*/ , "ParoiB" );


    M_exporter->step( time )->setMesh( __mesh );
    M_exporter->step( time )->add( "u", U.template element<0>() );
    M_exporter->step( time )->add( "p", U.template element<1>() );
    M_exporter->step( time )->add( "mesh_velocity", mesh_velocity );
    M_exporter->step( time )->add( "aux", aux );
    M_exporter->save();
#endif

    M_exporter->step( M_cpt_export )->setMesh( mesh_fluid );
    M_exporter->step( M_cpt_export )->add( "vitesse", U_fluid->template element<0>() );
    M_exporter->step( M_cpt_export )->add( "pression", U_fluid->template element<1>() );
    M_exporter->step( M_cpt_export )->add( "mesh_velocity", mesh_velocity );
    M_exporter->step( M_cpt_export )->add( "aux", aux );
    M_exporter->save();


    //M_exporter_ref->step( M_cpt_export )->setMesh( mesh_reference_fluid );
    M_exporter_ref->step( M_cpt_export )->setMesh( M_Ufluid_ref->mesh() );
    M_exporter_ref->step( M_cpt_export )->add( "u", M_Ufluid_ref->template element<0>() );
    M_exporter_ref->save();

    ++M_cpt_export;

    if (true) std::cout << "[NavierStokes] : Export NS finish\n";

} // NavierStokes::export


template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
void
NavierStokes<Dim, Order, Entity>::changeRepository(std::string pathLoc)
{

    this->application()->changeRepository( boost::format( "%1%/%2%/P%3%/h_%4%/"+ pathLoc )
                                           % this->application()->about().appName()
                                           % convex_type::name()
                                           % Order
                                           % this->meshSize
                                           );


}


template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
double
NavierStokes<Dim, Order, Entity>::diffL2( element_fluid_type& U1, element_fluid_type& U2 )
{
    using namespace vf;

  element_fluid_0_type u1 = U1.template element<0>();
  element_fluid_1_type p1 = U1.template element<1>();
  element_fluid_0_type u2 = U2.template element<0>();
  element_fluid_1_type p2 = U2.template element<1>();

  mesh_ptrtype mesh = U1.functionSpace()->mesh();

  double __norm = math::sqrt( integrate( elements(mesh),// _Q<2*(Order+OrderGeo-1)>(),
					 trans(idv(u1))*idv(u1) +idv(p1)*idv(p1) ).evaluate()(0,0) );
  double __diff = math::sqrt( integrate(elements(mesh),// _Q<2*(Order+OrderGeo-1)>(),
					trans(idv(u1)-idv(u2))*(idv(u1)-idv(u2))
					+(idv(p1)-idv(p2))*(idv(p1)-idv(p2)) ).evaluate()(0,0) );
  if ( __norm != 0 ) __diff /= __norm;

  return __diff;
}

} // end namespace NavierStokes

} // end namespace Feel

// instantiation
template class Feel::FSI::NavierStokes<2,2,Feel::Simplex>;
//template class Feel::NavierStokes<2,4,Feel::Simplex>;

//template class Feel::NavierStokes<3,2,Feel::Simplex>;

